GSSHA Workflow ¶
This notebook shows a sample workflow for running hydrology simulations using the GSSHA rectangularly gridded simulator, supported by a suite of primarily open-source Python libraries. The workflow consists of:
- Selecting parameters using widgets in a Jupyter notebook to control the model to simulate, including a watershed shape file.
- Visualizing the watershed shape in a geographic context (projected into a suitable coordinate system and overlaid on map tiles from a web tile server).
- If necessary, editing that watershed shape by hand and creating a new shape file with the edited result.
- Selecting parameters to control the simulation, potentially overriding some selected earlier for the model creation (e.g. if running numerous conditions as a parameter sweep).
- Visualizing and reviewing the inputs to the simulation.
- Running the underlying simulation, collecting data on flood depth at each time point as well as the overall maximum flood depth per grid cell.
- Visualizing the flood depth over time and the maximum flood depth.
- Analyzing the simulation speed to help shape expectations for computational requirements for future runs.
Each of these steps is configured directly in this notebook, and can thus easily be scripted or iterated as needed. The set of parameters and precisely how they are configured is still being improved, and it can likely be made into a better match to users' needs in this domain. This workflow relies on fast raster regridding added to Datashader and exposed via HoloViews as part of the EarthSim project.
The underlying environment needed to run this workflow is set up as described in the README , and though already functional will need to be greatly simplified to be more usable and maintainable in practice. The workflow currently relies on downloading data from external servers that can be slow to access from some parts of the internet, so you may see widely varying runtime speeds, especially the first time it is run.
from datetime import datetime, timedelta
import os
import glob
import shutil
import param
import paramnb
import numpy as np
import xarray as xr
import geoviews as gv
import holoviews as hv
import quest
import earthsim.gssha as esgssha
import earthsim.gssha.model as models
import cartopy.crs as ccrs
from earthsim.gssha import download_data, get_file_from_quest
from holoviews.streams import PolyEdit, BoxEdit, PointDraw, CDSStream
from holoviews.operation.datashader import regrid, shade
from earthsim.io import save_shapefile, open_gssha, get_ccrs
regrid.aggregator = 'max'
hv.extension('bokeh')
%output holomap='scrubber' fps=2
shutil.rmtree('./vicksburg_south',ignore_errors=True)
Configure model parameters ¶
model_creator = esgssha.CreateGSSHAModel(name='Vicksburg South Model Creator',
mask_shapefile='../data/vicksburg_watershed/watershed_boundary.shp',
grid_cell_size=90)
paramnb.Widgets(model_creator,initializer=paramnb.JSONInit())
Draw bounds to compute watershed ¶
Allows drawing a bounding box and adding points to serve as input to compute a watershed:
%%opts Polygons [width=900 height=500] (fill_alpha=0 line_color='black')
%%opts Points (size=10 color='red')
tiles = gv.WMTS('http://c.tile.openstreetmap.org/{Z}/{X}/{Y}.png',
crs=ccrs.PlateCarree(), extents=(-91, 32.2, -90.8, 32.4))
box_poly = hv.Polygons([])
points = hv.Points([])
box_stream = BoxEdit(source=box_poly)
point_stream = PointDraw(source=points)
tiles * box_poly * points
if box_stream.element:
element = gv.operation.project(box_stream.element, projection=ccrs.PlateCarree())
xs, ys = element.array().T
bounds = (xs[0], ys[1], xs[2], ys[0])
print("BOUNDS", bounds)
if point_stream.element:
projected = gv.operation.project(point_stream.element, projection=ccrs.PlateCarree())
print("COORDINATE:", projected.iloc[0]['x'][0], projected.iloc[0]['y'][0])
Inspect and edit shapefile ¶
The plot below allows editing the shapefile using a set of tools. The controls for editing are as follows:
- Double-clicking the polygon displays the vertices
- After double-clicking the point tool is selected and vertices can be dragged around
- By tapping on a vertex it can be selected, tapping in a new location while a single point is selected inserts a new vertex
- Multiple points can be selected by holding shift and then tapping or using the box_select tool
-
Once multiple vertices are selected they can be deleted by selecting the point editing tool and pressing
backspace
%%opts Shape [width=900 height=500 tools=['box_select']] (alpha=0.5)
mask_shape = gv.Shape.from_shapefile(model_creator.mask_shapefile).last
tiles = gv.WMTS('http://c.tile.openstreetmap.org/{Z}/{X}/{Y}.png')
vertex_stream = PolyEdit(source=mask_shape)
tiles * mask_shape
If any edits were made to the polygon in the plot above we save the
watershed_boundary.shp
back out and redisplay it to confirm our edits were applied correctly:
%%opts Shape [width=600 height=400] (alpha=0.5)
if vertex_stream.data:
edited_shape_fname = './vicksburg_watershed_edited/watershed_boundary.shp'
dir_name = os.path.dirname(edited_shape_fname)
if not os.path.isdir(dir_name): os.makedirs(dir_name)
save_shapefile(vertex_stream.data, edited_shape_fname, model_creator.mask_shapefile)
model_creator.mask_shapefile = edited_shape_fname
mask_shape = gv.Shape.from_shapefile(edited_shape_fname).last
mask_shape = mask_shape.opts() # Clear options
mask_shape
Configure simulation parameters ¶
sim = esgssha.Simulation(name='Vicksburg South Simulation', simulation_duration=60*60,
rain_duration=30*60, model_creator=model_creator)
paramnb.Widgets(sim,initializer=paramnb.JSONInit())
Create the model ¶
Note that the above code demonstrates how to collect user input, but it has not yet been connected to the remaining workflow, which uses code-based specification for the parameters.
if sim.model_creator.project_name not in quest.api.get_collections():
quest.api.new_collection(sim.model_creator.project_name)
paramnb.Widgets(sim.model_creator,initializer=paramnb.JSONInit())
# temporary workaround until workflow cleanup/parameterization is done
if sim.model_creator.project_name == 'test_philippines_small':
sim.model_creator.roughness = models.GriddedRoughnessTable(
land_use_grid=get_file_from_quest(sim.model_creator.project_name, sim.land_use_service, 'landuse', sim.model_creator.mask_shapefile),
land_use_to_roughness_table='./philippines_small/land_cover_glcf_modis.txt')
else:
sim.model_creator.roughness = models.GriddedRoughnessID(
land_use_grid=get_file_from_quest(sim.model_creator.project_name, sim.land_use_service, 'landuse', sim.model_creator.mask_shapefile),
land_use_grid_id=sim.land_use_grid_id)
sim.model_creator.elevation_grid_path = get_file_from_quest(sim.model_creator.project_name, sim.elevation_service, 'elevation', sim.model_creator.mask_shapefile)
model = sim.model_creator()
# add card for max depth
model.project_manager.setCard('FLOOD_GRID',
'{0}.fgd'.format(sim.model_creator.project_name),
add_quotes=True)
# Add time-based depth grids to simulation
"""
See: http://www.gsshawiki.com/Project_File:Output_Files_%E2%80%93_Required
Filename or folder to output MAP_TYPE maps of overland flow depth (m)
every MAP_FREQ minutes. If MAP_TYPE=0, then [value] is a folder name
and output files are called "value\depth.####.asc" **
"""
model.project_manager.setCard('DEPTH', '.', add_quotes=True)
model.project_manager.setCard('MAP_FREQ', '1')
# add event for simulation (optional)
"""
model.set_event(simulation_start=sim.simulation_start,
simulation_duration=timedelta(seconds=sim.simulation_duration),
rain_intensity=sim.rain_intensity,
rain_duration=timedelta(seconds=sim.rain_duration))
"""
# write to disk
model.write()
Review model inputs ¶
Load inputs to the simulation ¶
name = sim.model_creator.project_name
CRS = get_ccrs(os.path.join(name, name+'_prj.pro'))
roughness_arr = open_gssha(os.path.join(name,'roughness.idx'))
msk_arr = open_gssha(os.path.join(name, name+'.msk'))
ele_arr = open_gssha(os.path.join(name, name+'.ele'))
roughness = gv.Image(roughness_arr, crs=CRS, label='roughness.idx')
mask = gv.Image(msk_arr, crs=CRS, label='vicksburg_south.msk')
ele = gv.Image(ele_arr, crs=CRS, label='vicksburg_south.ele')
Shapefile vs. Mask ¶
tiles * regrid(mask) * mask_shape
Elevation ¶
tiles * regrid(ele) * mask_shape
Roughness ¶
tiles * regrid(roughness) * mask_shape
Run Simulation ¶
from gsshapy.modeling import GSSHAFramework
# TODO: how does the info here relate to that set earlier?
# TODO: understand comment below
# assuming notebook is run from examples folder
project_path = os.path.join(sim.model_creator.project_base_directory, sim.model_creator.project_name)
gr = GSSHAFramework("gssha",
project_path,
"{0}.prj".format(sim.model_creator.project_name),
gssha_simulation_start=sim.simulation_start,
gssha_simulation_duration=timedelta(seconds=sim.simulation_duration),
# load_simulation_datetime=True, # use this if already set datetime params in project file
)
# http://www.gsshawiki.com/Model_Construction:Defining_a_uniform_precipitation_event
gr.event_manager.add_uniform_precip_event(sim.rain_intensity,
timedelta(seconds=sim.rain_duration))
gssha_event_directory = gr.run()
Visualizing the outputs ¶
Load and visualize depths over time ¶
depth_nc = os.path.join(gssha_event_directory, 'depths.nc')
if not os.path.isfile(depth_nc):
# Load depth data files
depth_map = hv.HoloMap(kdims=['Minute'])
for fname in glob.glob(os.path.join(gssha_event_directory, 'depth.*.asc')):
depth_arr = open_gssha(fname)
minute = int(fname.split('.')[-2])
# NOTE: Due to precision issues not all empty cells match the NaN value properly, fix later
depth_arr.data[depth_arr.data==depth_arr.data[0,0]] = np.NaN
depth_map[minute] = hv.Image(depth_arr)
# Convert data to an xarray and save as NetCDF
arrays = []
for minute, img in depth_map.items():
ds = hv.Dataset(img)
arr = ds.data.z.assign_coords(minute=minute)
arrays.append(arr)
depths = xr.concat(arrays, 'minute')
depths.to_netcdf(depth_nc)
else:
depths = xr.open_dataset(depth_nc)
depth_ds = hv.Dataset(depths)
depth_ds.data
Now that we have a Dataset of depths we can convert it to a series of Images.
%%opts Image [width=600 height=400 logz=True xaxis=None yaxis=None] (cmap='viridis') Histogram {+framewise}
regrid(depth_ds.to(hv.Image, ['x', 'y'])).redim.range(z=(0, 0.04)).hist(bin_range=(0, 0.04))
We can also lay out the plots over time to allow for easier comparison.
%%opts Image [width=300 height=300 logz=True xaxis=None yaxis=None] (cmap='viridis')
regrid(depth_ds.select(minute=range(10, 70, 10)).to(hv.Image, ['x', 'y']).redim.range(z=(0, 0.04))).layout().cols(3)